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1 Introduction 



ABS methods were introduced by AbafFy, Broyden and Spedicato (1982/84) [|AbBS 82 



I AbBS 84 ] originally for solving systems of linear equations. The basic ABS class was later 
generalized to the so-called scaled ABS class and, subsequently, applied to linear-squares, 
nonlinear equations and optimization problems, see for instance, Abaffy and Spedicatio 



| AbSp 89 I, Spedicato [ Sped 97 |, [ Sped 99 |, Spedicato, Xia Z. and Zhang L. | SpXZ 00 



and Zhang L., Xia Z. and Feng E. |]ZhXF 99" |. In this paper we first review the general 
scheme for solving linear determined or undetermined systems. 

Let us consider the general linear systems, where rank(^) is arbitrary, 



or 



1, • • • ,m 



where 



^CORA, Department of Applied Mathematics, Dalian University of Technology Dalian 116024, China, 

xingli@student.dlut.edu.cn 

^Division of Software, Legend Computer Company, Beijing 100085, China, liuyingeQlegend. com.cn. 
^Department of Mathematics, University of Bergamo, Piazza Rosate 2, 24129 Bergamo, Italy 



1 



X G M^, b £ iR", m <n and A 



The steps of the unsealed ABS class of algorithms are defined as follows: 



Basic (Unsealed) ABS Class of Algorithms: HAbBS 84 | , Algorithml jAbSp 89 

pp. 21-22 

Begin Algorithm 

(A) Initialization. 

Give an arbitrary vector xi € iR", and an arbitrarily nonsingular matrix Hi € -R" 
Set i = 1 and iflag=0. 

(B) Computer two quantities. 
Compute 

Si — Hidi 

Ti = t'- ei = aj Xi - 



(C) Check the compatibility of the system of linear equations. 
If Si / then goto (D). 

If Sj = and = then set 

Hi+i = Hi 

and goto (F), the i-th equation is a linear combination of the previous equations. 
Otherwise stop, the system has no solution. 

(D) Computer the search vector pi £ iR" by 



Pi = Hjzi 



(1.1) 



where Zj, the parameter of Broyden, is arbitrary save that 



zjHitti / 



:i.2) 



(E) Update the approximation of the solution Xi by 



Xi+i = Xi- aipi 



:i.3) 
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where the stepsize Oi is computed by 

ai = Ti/aJpi (1.4) 

If i = m stop; a^m+i solves the system. 
(F) Update the (Abaffian) matrix Hi. Compute 

Hi+i = Hi- Hiaiwf Hi/wf Hiai (1.5) 

where Wi G IRP, the parameter of Abaffy, is arbitrary save for the condition 

wjHiai = 1 or 7^ (1.6) 



(G) Increment the index i by one and goto (B). 

We define n by i matrices Ai, Wi and Pi by 

Ai = {ai, ■■■ , aiY, Wi = {wi, ■ ■ ■ , Wi), Pi = {pi, ■ ■ ■ , Pi) (1-7) 
End Algorithm 

Some properties of the above recursion, see for instance, Abaffy and Spedicato (1989), 
| AbSp 89 I, are hsted below that are the basic formulae for use later on. 



a. Imphcit factorization property 

with Li nonsingular lower triangular. 

b. Null space characterizations 



AT Pi = Li 



M{H,^,) = n{Af), N{Hl,) = n{w,\ 



i+ll 
I+l) 



:i.9) 



where J\f= Null and 7^= Range, 
c. The linear variety containing all solutions to Ax = b consists of the vectors of the form 

X = xt+i + H^^^q (1.10) 

where q € iR" is arbitrary. 



Much progress on the computational aspect of ABS algorithms have been made since the 
ABS algorithms were found in the early 1980's. In the last several years lots of work has been 
done on enlarging, improving and completing ABSPACK that is a package of ABS algo- 



rithms in FORTRAN codes due to Spedicato and his collaborators |Bodo 89 ], |Bodo 90a 



i Bodo 00a|] , [podo 00b|| , jBodo 01a| ], i"BoLS 00a|] , [|BPLS 00a|] , [|BPLS 00b| ], jBPLS 01 



In the last two years we investigated possibility of establishing C-|--|-/VC-|--|- ABS software, 
i.e. ABS software written in C-|— I-/VC-I— 1-. Part of work in the subject is presented in this 
paper concerning the implicit LU algorithms. 

In this report a text of some codes implementing the implicit LU algorithms for solving linear 
determined and undetermined systems with n variables and m equations is given. The codes 
are written in C-|— |- language using Matrix class (with CMatrix as its class name) and Vector 
class (with CVector as its class name) that are made by ourselves for constructing the ABS 
software. The paper is organized based on the following three algorithms: 

1. Function iLUa: the implicit LU method, where A is regular (i.e. all principal subma- 
trices are nonsingular) , without pivoting for determined and undetermined systems. 

2. Function iLUaPivotC: the implicit LU method with column pivoting and explicit col- 
umn interchanges for linear determined and undetermined systems. Since the column 
pivot cause that the ordering of components of solution is changed, it is necessary to 
recover the ordering of components of solution and related codes are given. 

3. Function iLUaPivotR: the implicit LU method with row pivoting and explicit row 
interchanges for determined and undetermined systems. 

In the next section three schemes of implicit LU algorithms and some general properties are 
given. Codes implementing these algorithms are listed in section 3. Finally a main program 
to test algorithms is presented. 



2 Implicit LU Algorithms 



2.1 Schemes of implicit LU algorithms 



We give a short description of three versions of the implicit LU algorithms below: 
Implicit LU Algorithm Without Pivoting (iLUa): ||AbBS 82 || , ||AbBS 84 || , [^ped 01 
(iLUa is abbreviation for implicit LU algorithm without pivoting, under the assumption 



A 



of A is regular) 

Set xi = 0, Hi = I and i = 1 

For i = 1 to m do 

Set Si = Hitti 

di — 5^ 62 

Xj+i = Xi - {{afxi - bi)/di)Hjei 
if i < m, then set 

Hi+i = Hi — SiejHi/di 

enddo 



Implicit LU Algorithm With Column Pivoting (iLUaPivotC): jAbSp 89 | , jBodo 00a 
(iLUaPivotC is the abbreviation of "impUcit LU algorithm with column pivoting" and col- 
umn interchanges, without the regularity of A. The ordering of components of solution is 
changed because of the column pivot) 
Set xi = 0, Hi = I and i = 1 
For i = 1 to m do 
Set Si = HiUi 

(only (i — l)(n — i + 1) nonzero elements of Hi are used) 

Determine di = [sfe^J, such that \sfek^ \ = max{|sfej| \ j = i, ■ ■ ■ ,n} 

(only n — i + 1 nonzero elements of Sj are used) 

If ki ^ i, then swap columns of A and elements of Xj and Si with these indeices 
Set Xj+i = Xi- {{ajxi - bi)/di)Hfei 
(only i nonzero elements of Xi are updated) 
if i < m, then set -f^i+i = Hi — SicfHi/di 
(only i{n — i) nonzero elements of Hi^i are updated) 
enddo 

Implicit LU Algorithm With Row Pivoting (iLUaPivotR): 

(iLUaPivotR stands for implicit LU algorithm with row pivoting and explicit row inter- 
changes, without the regularity of A and change of components of solution vector) 
Set xi = 0, Hi = I and i = 1 
For i = 1 to m do 

Determine di = Ipfa^. |, such that bfflfcj = vaax{\pfaj \ \ j = i, ■ ■ ■ ,n} 
If ki ^ i, then swap rows of A and elements of Xi and b with these indices 
Set Si = HiUi 



(only (i — l)(n — i + 1) nonzero elements of are used) 
Set Xi+i = Xi- {{afxi - bi)/di)Hjei 
(only i nonzero elements of Xi are updated) 
ii i <m, then set -ffi+i = Hi — SieJ Hi/di 
(only i(ra — i) nonzero elements of -ffj+i are updated) 
enddo 



2.2 Properties of general implicit LU algorithms 

The general implicit LU algorithm is obtained by the parameter choices Hi = I, Zi = Si, 
Wi = Si- Some properties of this class of algorithms is listed as foUowings: 

(a) The algorithm is well defined iff A is regular (i.e., all principal submatrices are non- 
singular). Otherwise pivoting has to be performed. 

(b) Since W^Hi+i = [Jj, Oj-^iJj+i = 0, the first i rows of the Abaffian matrix must be 
zero. More precisely, the Abaffian matrix has the following structure, with Ki G 



H 



+1 







K T 



(c) Only the first i components of pi can be nonzero and the i-th is unity. Hence the 
matrix Pj is unit upper triangular, so that the implicit factorization A = LP~^ is of 
the LU type, with unit on the diagonal. 

(d) Only Ki has to be updated. The algorithm requires nm? — 2m^ multiplications plus 
lower-order terms. Hence, for m = n there are n^/3 multiplications plus low-order 
terms, which is the same cost as for the classical LU factorization or Gaussian elimi- 
nation (which are two essentially equivalent process). 



3 Codes of Implicit LU Algorithms 
3.1 ILUa code 

function BOOL iLUa(CMatrix m_A, CVector v_b, CVector &v_x, double epl, 

double ep2) 



// Under the condition of regularity of A, iLUa determines the solution 
// of the linear systems Ax = b ( A with dimension mxn, m<n), using 
// implicit LU algorithm without pivoting 

// the type of return value of function is BOOL, true if system has a 
// solution, otherwise the return value is false 

// m_A = an object of CMatrix class, denote coefficient matrix 

// v_b = an object of CVector class, denote right hand-side vector 

// v_x = an object of CVector class, denote solution vector 

// (using reference operator to output the solution vector) 

// epl = a double type value of dependency control parameter 

// ep2 = a double type value of residual control parameter 

{ 

try 
{ 

CVector v_xl; 

// declare an object of CVector class 
v_xl=v_x; 

CVector v_s; 

// declare an object of CVector class, denotes the vector Sj = i^jOj 
CMatrix m_H(l .0,m_A.m_cols) ; 

// declare an object of CMatrix class, denotes the Abaffian matrix, 
// the initial matrix is unit matrix 
double ns; 

// declare a double precision number, denotes , the norm of p 

double r; 

// declare a double precision number, ri = ajxi — bi 
int iflag=0; 

int i=l; 

// iteration 

while (i<=m_A . m_rows) 

{ 
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/ / compute Si = HiUi 
if(i==l) 

v_s=in_A.GetRow(i) .TransO ; 

else 

{ 

v_s[i-l]=0.0; 

f or (int il=i ; il<=m_A .m_cols ; il++) 
{ 

double suml=0.0; 
forCint jl=l; jl<i; 

suinl=suinl+m_H . Get Value (il , j 1) *m_A . Get Value (i , j 1) ; 
v_s [il] =suml+m_A . GetValue (i , il) ; 

} 



r=Dot(m_A.GetRow(i) , v_xl)-v_b [i] ; 
ns=v_s . Module ( ) ; 
if (ns<=epl) 
{ 

if (fabs(r)<=ep2) 

// i-fh. row is linearly dependent with first i — \ rows 
{ 

if lag=if lag+1 ; 
i=i+l; 

} 

else 
{ 

// system is incompatible 
iflag=-i; 

AfxMessageBoxC'No Solution"); 
break; 

} 

} 

else 
{ 
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// update solution Xi 
double temp=0; 
temp=r/v_s [i] ; 
for(int i2=l;i2<=i;i2++) 

v_xl [12] =v_xl [12] -temp*m_H . GetValue (1 , 12) ; 

// update projection matrix Hi 

If (l<m_A.m_cols) 
{ 

f or (Int 13=1+1 ; 13<=m_A.m_cols ; 13++) 
{ 

forClnt j3=l; j3<=l; j3++) 
{ 

double tt ; 

tt= (v_s [13] *m_H . GetValue ( 1 , j 3) ) /v_s [1] ; 
tt=in_H . GetValue (13 , j 3) -tt ; 
m_H . SetValue ( 13 , j 3 , tt ) ; 

> 

} 

for (13=1 ; 13<=1 ; 13++) 

m_H . SetValue (1 , 13 , . 0) ; 

} 

1++; 

} 

} 

If (lflag>=0) 
{ 

v_x=v_xl ; 
return (true) ; 

} 

else 

return(f alse) ; 

> 

catch(CErrorExceptlon *e) 
{ 



Q 



Af xMessageBox(e->GetErrorInf o ) ; 
e->Delete() ; 
exit(l) ; 

} 

} 

3.2 iLUaPivotC code 

function BOOL iLUaPivotC (CMatrix m_A, CVector v_b, CVector &v_x, double epl, 

double ep2) 

// iLUPivotC determines the solution of the linear systems Ax = b (A with 
// dimension mxn, m < n ) using implicit LU algorithm with column pivoting 
// and explicit column interchanges, and the order of components solution 
// vector are changed 

// the type of return value of function is BOOL, true if system has a solution 
// otherwise the return value is false 

// m_A = an object of CMatrix class, denote coefficient matrix 

// v_b = an object of CVector class, denote right hand-side vector 

// v_x = an object of CVector class, denote solution vector 

// (using reference operator to output the solution vector) 

// epl = double type value of dependency control parameter 

// ep2 = double type value of residual control parameter 

{ 

try 
{ 

CVector v_xl; 

// declare an object of CVector class 
v_xl=v_x; 

CVector v_s; 

// declare an object of CVector class, denotes the vector Si = HiUi 
CMatrix m_H(l .0,m_A.m_cols) ; 

// declare an object of CMatrix class, denotes the Abaffian matrix, 

in 



// the initial matrix is unit matrix 
int ki ; 

// declare an integer number, denotes the ordering of pivot 
double r; 
int iflag=0; 
double d; 

int* Index; 

// declare a pointer to type integer 
Index=new int [m_A.m_cols] ; 

// use new operator to allocate a pointer to an array of integer with 

// dimension m_A.m_cols 

for (int ii=l ; ii<=m_A.m_cols ; ii++) 

Index [ii] =ii ; 
// initialize the array of integer 
int i=l; 

// iteration 

while (i<=m_A.m_rows) 

{ 

// compute Si = HiUi 

if(i==l) 

v_s=m_A.GetRow(i) .TransO ; 
else 
{ 

v_s[i-l]=0.0; 

f or(int il=i; il<=m_A.m_cols; il++) 
{ 

double suml=0.0; 
for (int jl=l; jl<i; 

suinl=suml+m_H . Get Value (il , j 1) *m_A . Get Value (i , j 1) ; 

v_s[il]=suinl+m_A. Get Value (i,il) ; 

} 

} 

1 1 



// pivoting 
d=-l; 

f or (int j 2=i ; j 2<=m_A . m_cols ; j 2++) 
{ 

if (d<fabs(v_s[j2])) 
{ 

d=fabs(v_s[j2]) ; 
ki=j2; 

} 

} 

if (ki!=i) 

// swap i-th. column and fc^-th column of A 

// swap i-th component and ki-th component of Xj and Sj 

// save ordering of pivot in Index 

{ 

m_A=m_A . SwapCol (ki , i) ; 

v_x=v_x . SwapElement (ki , i) ; 

v_s=v_s . SwapElement (ki , i) ; 

int it ; 

it=Index [i] ; 

Index [i] =Index [ki] ; 

Index [ki] =it ; 

} 

r=Dot(m_A.GetRow(i) ,v_xl)-v_b[i] ; 

if (f abs (v_s [i] ) <=epl) 

{ 

if (fabs(r)<=ep2) 

// i-th row is linearly dependent with first i — 1 rows 
{ 

if lag=if lag+1 ; 
i=i+l; 

} 

else 
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// system is incompatible 
{ 

iflag=-i; 

AfxMessageBoxC'No Solution"); 
break; 

} 

} 

else 
{ 

// update solution Xi 
double temp=0; 
temp=r/v_s [i] ; 
for(int i2=l;i2<=i;i2++) 

v_xl [12] =v_xl [12] -temp*m_H . GetValue (i , 12) ; 

// update projection matrix 

if (i<m_A.m_cols) 

{ 

f or (int 13=1+1 ; i3<=m_A .m_cols ; 13++) 
{ 

for (int j3=l; j3<=i; j3++) 
{ 

double tt; 

tt= (v_s [13] *m_H . GetValue (1 , j 3) ) /v_s [1] ; 
tt=m_H . GetValue (13 , j 3) -tt ; 



m_H.SetValue(i3, j3,tt) ; 

} 

} 

for (13=1 ; 13<=1 ; 13++) 
m_H . SetValue (1 , 13 , . 0) ; 



> 

1++; 

> 

if (lflag>=0) 
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// exchange components of solution to adapt the original system 
{ 

f or (int i4=l ; i4<=m_A .m_cols ; i4++) 
{ 

for (int j4=l; j4<=m_A.m_cols; j4++) 
{ 

if (Index [j 4] ==i4) 
v_x[i4]=v_xl[j4] ; 

} 

} 

return (true) ; 

} 

else 

return(f alse) ; 

} 

catch(CErrorException *e) 
{ 

AfxMessageBox(e->GetErrorInf o() ) ; 
e->Delete() ; 
exit (1) ; 

} 

} 

3.3 iLUaPivotR code 

function BOOL iLUaPivotR (CMatrix &m_A, CVector &v_b, CVector &v_x, double epl, 

double ep2) 

// iLUPivotC determines the solution of the linear systems Ax = b (A with 
// dimension mxn, m <n ) using implicit LU algorithm with column pivoting 
// and explicit column interchanges 

// the type of return value of function is BOOLtrue if system has a solution 
// otherwise the return value is false 

// m_A = an object of CMatrix class, denote coefficient matrix 
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// v_b = an object of CVector class, denote right hand-side vector 

// v_x = an object of CVector class, denote solution vector (using reference 

// operator to output the solution vector) 

// epl = double type value of dependency control parameter 

// ep2 = double type value of residual control parameter 



try 
{ 

CVector v_xl ; 

// declare an object of CVector class 
v_xl=v_x; 



CVector v_s; 

// declare an object of CVector class, denotes the vector Si = HiUi 
double r; 

// declare a double precision number, ri = afxi — bi 
CMatrix m_H(l .0,m_A.m_cols) ; 

// declare an object of CMatrix class, denotes the Abaffian matrix, the 
// initial matrix is unit matrix 
int ki; 

// declare an integer number, denotes the ordering of pivot 
double d; 
int iflag=0; 
int i=l; 



// iteration 

while (i<=m_A.m_rows) 

{ 

// pivoting 
double mpt=0 
d=-l; 

f or ( int j 2=i ; j 2<=m_A . m_rows ; j 2++) 

1 



{ 

mpt=Dot(ni_H.GetRow(i) ,m_A.GetRow(j2)) ; 
if (d<fabs(mpt)) 

{ 

d=f abs (mpt) ; 
ki=j2; 

} 

} 

if (ki!=i) 

// swap z-th row and /cj-th row of A 

II swap z-th component and /cj-th component of Xi and Sj 

// save ordering of pivot in Index 

{ 

m_A=m_A.SwapRow(ki,i) ; 
v_x=v_x . SwapElement (ki , i) ; 
v_b=v_b . SwapElement (ki , i) ; 

} 

mpt=Dot(m_H.GetRow(i) ,m_A.GetRow(i)) ; 

// compute Si = Hiai 
if(i==l) 

v_s=m_A.GetRow(i) . Trans () ; 
else 
{ 

v_s[i-l]=0.0; 

f or(int il=i; il<=m_A.m_cols; il++) 
{ 

double suml=0.0; 
for(int jl=l; jl<i; jl++) 

suml=suml+m_H . GetValue (il , j 1) *m_A . GetValue (i , j 1) ; 
v_s[il]=suml+m_A.GetValue(i,il) ; 

} 

} 

r=Dot(m_A.GetRow(i) ,v_xl)-v_b[i] ; 



if (fabs(nipt)<=epl) 
{ 

if (fabs(r)<=ep2) 

// i-th. row is linearly dependent with first i — 1 rows 
{ 

if lag=if lag+1 ; 
i=i+l; 

} 

else 

// system is incompatible 
{ 

iflag=-i; 

Af xMessageBoxC'No Solution"); 
breaJs; 

} 

} 

else 
{ 

// update solution 
double temp=0; 
temp=r/mpt ; 

for(int i2=l;i2<=i;i2++) 

v_xl [i2] =v_xl [i2] -temp*m_H . Get Value (i , i2) ; 
// update projection matrix Hi 
if (i<m_A.m_cols) 
{ 

f or (int i3=i+l ; i3<=m_A .m_cols ; i3++) 
{ 

for (int j3=l; j3<=i; j3++) 
{ 

double tt; 

tt= (v_s [i3] *m_H . Get Value (i , j 3) ) /mpt ; 
tt=m_H . Get Value (i3 , j3) -tt ; 
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m_H.SetValue(i3, j3,tt) ; 

} 

} 

for(i3=l;i3<=i;i3++) 

m_H. Set Valued, 13, 0.0) ; 

} 

1++; 

} 

} 

If (lflag>=0) 
{ 

v_x=v_xl ; 
return (true) ; 

} 

else 

return(f alse) ; 

} 

catch(CErrorException *e) 
{ 

AfxMessageBox(e->GetErrorInf o() ) ; 
e->Delete() ; 
exlt(l) ; 

} 

} 

4 A Main Program for LU Algorithm 

To make ABS algorithms be used more widely, we encapsulate Matrix class, 
Vector class and the ABS algorithms modules into the library modules 
ABSDLL.dll. A main program for the use of module(function) iLUaPivotC 
that solves linear system Ax = b (with Micchelli-Fiedler matrix as coefficient 
matrix) is as following: 

#lnclude "stdafx.h" 

IS 



#include "absalg.h" 
#include "iostream.h" 



void mainO 
{ 

try 
{ 

int m=1000; 
int n=1000; 

HINSTANCE hLib=AfxLoadLibrary(" ABSDLL.dll") ; 

CMatrix a(m,n) ; 
CVector b(m) ; 

for (int i=l ; i<=m; i++) 
{ 

for(int j = l; j<=ii; 
{ 

a.SetValue(i, j ,abs(i-j)) ; 

} 

> 

for (int k=l;k<=in;k++) 
b[k]=k; 

CVector xx(0.0,n) ; 

if ( ! iLUaPivotC (a , b , xx , 1 . Oe-7 , 1 . Oe-7) ) 
{ 

cout«("Do you want a least-square solution?\endl") ; 

char c; 

cin»c; 

if (c=='y' I |c=='yO 

// calling the modules for solving the least-squares problem 
else 



IQ 



exit(O) ; 

} 

else 
{ 

CVector r; 
r=a*xx-b ; 
er=r. Module ; 

cout«("errorbound=°/of \n" ,er) ; 
Af xFreeLibrary(hLib) ; 

} 

} 

catch(CErrorException *e) 
{ 

Af xMessageBox(e->GetErrorInf o() ) ; 
e->Delete() ; 

} 

} 

Remarks 

As for numerical experiments, we liave made some tests, for example, matrices with elements 
dij = {i — j\, 1 < i < m, 1 < j < n (Micchelli-Fiedler matrix) and matrices with elements 
aij = \i — jp, l<i<m, l<j<n. The results show that the algorithms are efficient. 
Further numerical experiments are in progress. 
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